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ABSTRACT 


This report documents a comparative study of the 
problems encountered in implementing the kinematic 
equations of rigid-body rotation. Comparing the efficiency 
and accuracy of the Euler angles and quaternions, the study 
indicates that quaternions are far superior to Euler angles 
for a large-angle simulation. The report explores the 
mechanization of the quaternion-constrained equations in 
order to obtain improved accuracy and simulation speed. It 
also describes the development of a new constraining tech- 
nique (called derivative constraint) for the quaternions and 
stability criteria for stable solution of the constrained 
equations. By use of these stability criteria, the optimum 
feedback gain constant for the constrained kinematic equa- 
tions can be selected. 
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DIGITAL SIMULATION OF ROTATIONAL KINEMATICS 


by 

Ai Chun Fang and Benjamin G. Zimmerman 
Goddard Space Flight Center 


INTRODUCTION 

Because of the increased need for dynamic simulation of the rotational motion of rigid bodies, 
digital simulations have become a necessary process. To describe the angular orientation of a 
rigid body in such simulations, the use of both Euler angles and quaternions is very common. 
Therefore, it is worthwhile to make a comparative simulation study of mathematical models, using 
Euler angles and quaternions for the kinematic representations. 

This report presents the results obtained from such a study. The report includes: 

1. A discussion of Euler angles and quaternions as kinematic representations. A review of 
these techniques is given specifying their advantages and shortcomings. 

2. An examination of integration techniques. The fourth order Runge Kutta fixed step inte- 
gration technique and the variable step Runge Kutta Merson integration technique are selected for 
this study. Simplicity and efficiency are the primary considerations used in evaluating the relative 
merits of each method. 

3. Presentation of significant parameters for the two representations. The following three 
parameters are studied in this paper: 

(a) Simulation speed 

(b) Attitude error (a measure of the deviation between the computed solution and the analyt- 
ically determined solution) 

(c) Constraint error (a measure of how the computed quaternion transformation matrix main- 
tains its orthogonality) 

Although the first of these parameters depends on the computing system used, a relative com- 
parison can be made from the result presented,, Double precision check runs seem to indicate that 
the attitude and constraint errors measured result primarily from the formulation of the model 
and selection of integration technique and thus are not related to the computing system. 

4. Discussion of constraint error reduction techniques. In general speed and accuracy satisfy 
a reciprocal relationship; that is, an increase in the simulation speed means a reduction in the 
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accuracy. Since many simulations involve solutions extending over large time intervals, a tech- 
nique for restricting error growth without reducing simulation speed is extremely desirable. Two 
known techniques for reducing the constraint errors are reviewed and investigated and a new tech- 
nique is suggested and discussed. 

5. Determination of the stability condition. The stability of the constraint error reduction 
technique depends upon the selection of the gain constants employed. The empirical determination 
of the stability conditions for a given technique and the representation of these conditions as a 
simple function of the system T s parameters require extensive review of the computer solutions. 
Relationships are developed for both of the techniques for which stability is a consideration. 

6. Theoretical verification for the stability conditions. The stability conditions first deter- 
mined by observation of the computer solutions are then verified by the development of theoretical 
proofs. The conditions are thus mathematically justified for insuring stability of the computed 
solution and are easily implemented on a digital computer. 

BENCHMARK PROBLEM 

The mathematical model used in this study will be referred to as the benchmark rotational 
kinematics problem. The selection of the benchmark problem and the techniques used resulted 
from a trade-off between simplicity and applicability. 

The benchmark problem, illustrated in Figure 1, consists of a fixed, right-handed orthogonal 
axis system 1 F , 2 F , 3 F and a moving axis system 1 M , 2 M , 3 M , whose initial orientation 1 M0 , 2 M0 , 
3 M0 is obtained from l F , 2 F , 3 F by a rotation about 2 F through an angle B o . The moving axis 
system moves with a constant rate a> about axis 3 F . Thus 

Ip * l«o r ( sinB o)2 F . 

‘ = -3 f . (1) 

The actual orientation of the coordinate vectors of the moving system at time t with respect 
to the fixed coordinate system is 

1„ - (cos cot cos B , sin wt cos , “ sin B ) , 

M \ o o o / 

2 m “ (- sina;t, cos art, 0) , 

3 , r (cos '4 sin B , sin sin B , cos B ) . (2) 

M \ o o o I 

The attitude error is taken to be the maximum of the angular errors between the computed position 
of the coordinate axes and their actual position; that is, 
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e T sin 1 -|max JJ 1 M (actual) x 1 M (computed) | , 

| 2 m ( actual ) x 2 M ( computed) | , 

(3) 

In determining simulation speed, the programs were rerun with the error-computing portions 
removed. All programs were written in* FORTRAN; however, "good programming practices 1 ’ were 
used to minimize computer time (e.g., avoidance of unnecessary subroutines, avoidance of sub- 
scripted variables, etc.). The simulation speed presented is calculated from solutions generated on 
a SDS 9300 computer. 


3 M (actual) 


3 m ( computed) 
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Except where noted, the equations are considered to be normalized so that <*> - 1 cycle per sec. 
Thus both angles and the independent variable time are measured in terms of cycles (1 cycle = 2 ? t 
radians). Computer time is expressed in seconds and simulation speed is determined by dividing 
cycles of the solution by the actual computer time required to generate it. 

EULER ANGLE KINEMATICS 

The use of sequences of Euler angles as kinematic parameters is quite common since they are 
easy to visualize and are a direct mathematical model of physical gimbal mounts. They have dis- 
advantages in that both the kinematic differential equations and the transformation matrix involve 
sines and cosines of the Euler angles, the generation of which is time consuming. Also the kine- 
matic differential equations have a singularity for certain values of the Euler angle (i.e., for those 
orientations for which "gimbal lock” occurs in physical gimbals). The particular Euler angle se- 
quence selected for study is of the successive type; however, the results are general since any 
Euler angle sequence can be obtained from any other sequence by a suitable renumbering of the 
axis system, accompanied by suitable changes in the positive sense of the axes to maintain right - 
handedness 0 For Euler angles A, B, C about axes 1, 2, 3 (Figure 2) the kinematic differential equa- 
tions are 



Figure 2— Euler angles. 
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The singularity in this system occurs when B = ± u/2. The transformation matrix giving the moving 
system with respect to the fixed system is 


cos B cos C sin A sin B cos C + cos A sin C - cos A sin B cos C + sin A sin C 

-cosBsinC - sin A sin B sin C + cos A cos C cos A sin B sin C + sin A cos C 
sinB -sin A cos B cosAcosB 


(5) 


Evidently, the terms of this matrix correspond to Equation 2, and thus 

A - tan -1 (- sin wt tan B o ) , 

C = tan -1 (tan oJt/cos B o ) . 

B sin" 1 (cosa;t sin Bj . (6) 

As cot increases, B varies between +B q and -B o . Thus the value of B o determines how 
closely the solution approaches the singular value. 


FIXED-STEP INTEGRATION 

Fixed-step integration techniques have the advantage of simplicity and a possible speed 
advantage for problems of roughly constant harmonic content. For small B q values, A , B, C are 
approximately sinusoidal, making fixed-step integration practical. When B o is near ± n/ 2 , A and 
C approach square waves, requiring an extremely small integration step to maintain accuracy. 

Figure 3 shows attitude error, using the Runge Kutta technique versus integration step size 
and critical angle B o . This error represents the difference between a digital simulation of the 
benchmark problem and its true analytical solution. The estimated error per step for the inte- 
gration technique is of the order of h 4 , which gives an error per cycle of the order h 3 ; the curves 
of constant B have this form. 

O 

The selection of integration step size is accomplished best by trial and error; however, when 
a prior selection is required, the use of some fraction (e.g., 1/10 or less) of the period associated 
with the highest frequency to be encountered gives a reasonable estimate. 


VARIABLE-STEP INTEGRATION 

Variable-step integration techniques utilize an expression for the estimated error per step as 
a criterion for the dynamic adjustment of the step size. Conceptually, when the variation of the 
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solution is relatively smooth, the step size is increased, giving an increase in simulation speeds; 
when the variation of the solution increases, the step size is reduced, causing a reduction in simu- 
lation speed but tending to preserve the accuracy of the solution,, The variable- step integration 
technique used here is the Runge Kutta Merson technique (Reference 1), which has an error per 
step of the order of h 4 as in the previous technique but which requires five evaluations of the 
derivatives per step compared to four for the previous technique. Thus, for the same step size h, 
the variable-step technique will take approximately 25 percent longer computing derivatives, in 
addition to the time required for computing and evaluating the error and adjusting the step size. 

In order to adjust the step size to control the error, the following procedure is used. An esti- 
mated error for each of the state variables is computed, and the system estimated error e s is 
arbitrarily taken to be the maximum of the state variable estimated errors. If e s is greater than 
some preassigned value e max , then the solution moves back one step and proceeds, using a step 
size one-half as large. If e s is less than some preassigned value e min , then the step size is 
doubled in succeeding steps. If e min < e s < e max , then the solution proceeds with no change in 
step size. In order to keep the number of parameters to a minimum, e min was set equal to e max . 

If e max is greater than e min , then some increase in simulation speed can be expected; however, 
in no case will it be greater than 2. 

Although variation of the integration error criterion can produce sizable jump discontinuities 
in simulation speed when using variable-step integration, if a sufficient density of solutions is 
obtained, the data can be smoothed so as to be presentable in a form analogous to Figure 3. 
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Figure 4 presents the smoothed results for the benchmark problem using the Runge Kutta 
Merson integration technique. Attitude error rate is plotted versus critical angle B o and integra- 
tion speed, with the step size criterion indicated parametrically. Since the step size criterion is 
applied to the Euler angles rather than to the attitude error, some variation in the error criterion 
curves may be expected; however, they do indicate the predicted trend; i.e., as the critical angle 
increases, the simulation speed decreases, while the attitude error rate remains roughly constant. 
A comparison of Figure 3 and Figure 4 indicates that a variable integration technique has no sig- 
nificant advantage in simulations where the critical angle is a problem. 



Figure 4— Euler angles with variable-step integration. 


7 



QUATERNIONS 


When solving the dynamic problems associated with a rotating rigid body, the use of quaternions 
(or Euler's symmetric parameters) as kinematic parameters has the capacity of handling unre- 
stricted rotation since the singularity is eliminated completely in this system. Thus the quaternions 
have the proven advantage of allowing simulation of a tumbling body. Also since the equations do 
not involve trigonometric functions, an increase in simulation speed is indicated. However, since 
four parameters are computed rather than the three required to describe the attitude of a body, 
the techniques for handling the constraint equation that relates the parameters bears consideration. 
For the selected benchmark problem mentioned earlier, the desired rotation is accomplished by 
two rotations— one through the angle B o about the 2 F axis, and one through the angle &t about the 
axis 3 f . In quaternion notation, the first rotation, the second rotation, and the resulting rotation 
may be presented by q B , q c , and q, respectively. The quaternion q is represented 

q r e o + e ! ’ e 2^2 4 e 3^3 ■ (7) 

where e Q , e l , e 2 , e 3 are real numbers regarded as the coefficient of the basis quaternions (1, q L , 
q 2 , and q 3 ). These basis quaternions obey the following rules of multiplication: 


V ^ q 2 2 -- q 3 2 = " 1 


q i q 2 

" q 2 C ll 

= q 3 ' 



' q 3 q 2 

= q i ’ 


^3 *1 

" q i q 3 

= q 2 ' 

(8) 


If the attitude of the coordinate system being described is considered as being obtained from 
the reference coordinate system by a single rotation, then the quaternion can be intrepreted geo- 
metrically as consisting of a vector part with components e , e , e 3 having direction associated 
with that of the rotation axis and having magnitude equal to the sine of one-half the angle of rota- 
tion and scalar part e Q equal to the cosine of one-half the rotation angle. Thus the quaternions q B 
and q c expressed in the fixed -coordinate system are 


q B - cos B o /2 + q 2 sinB o /2 


and 


q c - cos cot / 2 + q 3 sina)t/2 


(9) 
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From the algebra of quaternions, q is then given by 


q " <lp Qp 


(cos cot/2 + q 3 sin cvt/2) (cos B o /2 + q 2 sin B q /2) 


Therefore, 


e Q - cos cot/2 cos B o /2 


s in cot/2 s in B /2 » 


e - cosojt/2sinB /2 , 


e 3 - sin cot/2 cos B q /2 • 


In terms of e Q , e p e 2 , and e 3 , Equation 2 is equivalent to 


[e o 2 t e 2 - q 2 — q 2 ' 2 ( Bl e 2 +e 0 c 3 ) , 2(e t ? 3 ~ e fl e 2 )] 

[ 2 ( C 1 e 2" e 0 C 3)’ C 0 +e 2 2_e i 2 _e 3 2 ' 2 ( e 2 e 3 + % C l)] 


[ 2 ( e i e 3 + e 0 e 2 )' 2 ( e 2 e 3 _e 0 e i)' e 0 + 


e 2 “ e 2 - e 2 
3 12 


] 


(10) 


( 11 ) 


( 12 ) 


The quaternion components are not all independent, but they satisfy the constraint equation 


e2 + p2t02 + 02 | 

0 1 2 3 


(13) 


Differentiating Equation 13 with respect to time gives 


e o e 0 f e i e ! + e 2 e 2 + e 3 e 3 = 0 


(14) 


The angular velocity components in the moving system can be represented as 


K ) = 2 ( e 0 *1 - e i % _e 3 ®2 + 6 2 4 3 ) ’ 

("mJ = 2 ( e 0®2 _e 2 4 0 “ e i 4 3 + e 3 4 l) ' 

K) = 2 ( 6 0 ®3 _ e 3% _e 2 «1 +e i * 2 ) • 


(15) 
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Solving these equations and Equation 14 simultaneously yields 


*0 = ( e i° J M l +e 2 “M 3 +e 3“M 3 ) • 

*1 = T ( e 0 W M 1 ' e 3\ +e 2“M 3 ) > 

e„ = “T fe. + e„ “ e, \ , 

2 2y3Mj 0 M 2 1 M 3 J 

^3 = -T( e 2 - Ml - e i-M 2 - e 0-M 3 ) ' (16) 

The numerical solution of these four differential equations can be obtained without difficulty 
from the digital computer. A fixed-step four- point Runge Kutta integration technique was used in 
writing a program for this purpose. The initial values of the quaternion components at time t(0) = 0 
are simply 


e o = cos ( B o/2) - 

e i = 0 , 

e 2 = Sin ( B </ 2 ) ' 

e 3 = 0 • (17) 

With any numerical integration technique, a computer gives only an approximation of the ideal 
solution. As noted before, machine-related errors (i.e 0 , finite word length) and software errors 
(i 0 e 0 , truncation of series representations of functions) are small compared with the error introduced 
by the integration technique. The integration step size h contributes errors that depend on the 
magnitude of h and tend to accumulate as the solution progresses. 

For convenience, errors are divided into two types, attitude error and constraint error. The 
latter is defined as 


€ 


(18) 


For the benchmark problem, both types of errors are independent of the angle B o . They are 
functions of step size h only. If h ° t - t # is measured in cycles, then the constraint error rate, 
an increment of e per cycle, is expressed 


Ae 



(19) 
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The constraint error will accumulate in time to an intolerable value. In one sense, a constraint 
error is more serious than a position error since it represents nonorthonormality in the associated 
transformation matrix— a physically unrealizable situation. Reducing the constraint error by re- 
ducing h requires a correspondingly larger number of integration steps and a longer computing 
time. Other techniques for reducing the accumulation of constraint error are certainly needed, 
especially when the calculations involve direct 
computation of spacecraft positions over a 
long period of time in which a small h is 
impractical. 

The plots of constraint error rate and at- 
titude error rate versus h are given in Fig- 
ure 5. The straight line shows the simulation 
time per cycle of cot required for various val- 
ues of h using the SDS 9300 computer. 

Three methods of reducing the constraint 
error are discussed in the following para- 
graphs; they are algebraic constraints, normalized constraints, and derivatives constraints. These 
methods can improve the accuracy of solutions over the range of h for which the original system 
is stable. 


CONSTRAINT ERROR RATE (cycle" 1 ) 

ATTITUDE ERROR RATE (cycle/cycle) 


io 8 io 7 io 6 io 5 io 4 



Figure 5— Quaternions with no constraint. 


ALGEBRAIC CONSTRAINT 


In recent years, this method has been practiced in analog simulations of spacecraft motion 
and has been instrumental in correcting the constraint errors (Reference 2). The method consists 
of adding to each equation in Equation 16 a term proportional to £e 0 , ee p ee 2 , and ee 3 , respec- 
tively, giving the following constrained equations: 


e n ~ ~ \ l e , + e f e, ^ + Kee_ » 

0 illMj 2 M 2 3M 3 j 0 

*1 = {( e o"M,- e 3\ +e 2\) +K£e i ’ 

4 2 " T ( e 3"M, + e 0%- e i" M3 ) + K€e 2 * 

4 3 = ( e 2 a, M 1 " e i"M 2 - e 0 &, M 3 ) + Kee 3 ' 


( 20 ) 


When Equation 20 is solved on a digital computer, the effectiveness of the technique for a given 
value of h requires a judicious choice for the constant K. A Liapunov function exists for the system 
represented by Equation 20, and, therefore, it is always stable with all k > 0. It can be shown that 
the constraint error is minimized when K is increased. However, because of quantization and 
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errors introduced by the integration technique , there is an upper bound of K above which Equa- 
tion 20 would become unstable. 

Until now, no criterion has existed for finding the maximum value of k where the constrained 
Equation 20 would still be stable for a given integration step size h (Reference 3). Such a 
criterion is derived in the following paragraphs. 

The relationship between h and optimum K for the benchmark problem was first obtained em- 
pirically by repeated computer solution of Equation 20 using fixed- step Runge Kutta integration. 
Since this integration technique requires evaluation of the equations four times for each step, it 
was decided that Equation 18, required to obtain e , would not be included in the integration loop, 
but only evaluated at the beginning of each step. It was observed from the solutions that for all 

values of K > 0 for which stable solutions were 
obtained, the constraint error did not increase 
with time, but assumed a constant value. As K 
was increased, this constant constraint error 
decreased, and for each h there was an optimum 
value K o of k, above which instability occurred. 
These results are presented in Figure 6. Com- 
parison with Figure 5 reveals that the initial 
rate of increase of the attitude error is not af- 
fected by the improvement of the constraint 
performance. However, since the constraint 
error remains constant, the technique will not 
experience the errors introduced by nonortho- 
normality as would become apparent in the pre- 
viously described approach as time increased. 

The upper bound of K where the constrained Equation 20 remained stable was found to be 


CONSTRAINT ERROR 

-ATTITUDE ERROR RATE ( cycle/cycle ) 



0 10 20 30 40 50 60 70 80 90 100 


VALUE OF K ABOVE WHICH TECHNIQUE IS UNSTABLE 

Figure 6— Quaternions with algebraic constraint. 


hK < 1 . (21) 

Additional solutions for non-normalized values of w revealed that the relation is independent 
of co and that K o depends on step size h only. The important fact that K o is not a function of angular 
speed oj is quite striking. Generally speaking Equation 21 is true for any kinematic system for 
which quaternions are used. This result can be inferred theoretically in the following way. 

If the first equation of Equation 20 is multiplied by e Q , the second equation by e l? the third 
equation by e 2 , the fourth equation by e 3 , and these equations are summed, then the terms in the 
righthand side that contain co M , and cancel each other. Consequently, we have 

e 0 ®0 + e i®l + e 2 ®2 + e 3 4 3 = K£ ( e 0 2 +e i 2+ e 2 2 +e 3 2 ) • (22) 
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Since 


+ e 2 + e 2 + 


(23) 


and 


e + e, e, + e + e„ e o - - k / 2 , 

00 11 22 33 


(24) 


Equation 22 can be written as 


-e/2 = Ke(l-e) 


or 


- k £ 2KC 


(25) 


By using the relation 


( 6 ) t +h ( € ) t _ A£ 

h ' h ’ 


Equation 25 becomes 


-Ac £ 2 Khe . (26) 

For a discretely sampled system to be stable, it is necessary that |e | be decreasing for in- 
finitely many integration steps; that is, 


l< £ >. 


k e >t 


which implies that Ae is of opposite sign to c and 

l A£ l = |(O t+h -( e >,| * |( £ ) t+h | + |<O t | 1 z|(e) t | • . (27) 

If Equation 27 is applied to Equation 26, Equation 21 results. 

This condition may require that K cannot be a very large number in some cases. For instance, 
if co is small, then a large value of step size h is indicated, and the maximum permissible K will be 
small. Hence K does not depend on o> directly but is affected indirectly by « through its effect on h. 
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NORMALIZED CONSTRAINT 


When the constraint error is present, the sum of the square of the components of a quaternion 



• is no longer equal to unity. The quaternion can be renormalized by dividing each of its components 
by this norm. This has suggested the idea that the normalization of the quaternion at each step of 
the integration might reduce the constraint error. 

The components of the renormalized quaternion have the forms 


2 + e 2 + e 2^ 1/2 ^ e 0 (l 4 l/ 2£ ) 

e 2 e 3 


(e 2 +e 2 +e 2 + e 2 V 

Vo 1 2 3 ) 


2 + e 2 + e 2^ 1/ 2 ^ e i( l + l/ 2£ ) 


fe 2 + e 2 + e 2 + e 2 V 

Vo 1 2 3 I 


fe 2 + e ■ 
Vo 1 


2 + e 2 + e 2 y/2 


e (1 + 1/2 6 ) 


and 


( e 0 2+ e i 2 + e 2 2 + e 3 2 ) 1/2 


e, (1 + 1/2 e) 


(28) 


Thus 1/2 <=e 0 , 1/2 ee,, 1/2 ee 2 , and 1/2 ee 3 may be considered as the additional terms being 
added to e 0 , e x , e 2 , and e 3 , respectively. The derivatives are then calculated when e 0 , e 1 , e 2 , and 
e 3 are replaced by e Q +1/2 fe 0) ej + 1/2 eej , e 2 + 1/2 ee 2 , and e 3 + 1/2 ee 3 in Equation 16. 

Accordingly, the quaternion is normalized once every time interval h to compute the approxi- 
mation of its components that are substituted and integrated in Equation 16. It is hoped that in this 
way the constraint error may be diminished in each interval, and the accuracy of the final so- 
lution can be improved. 

The result of the computation indicated that the situation is not that simple. On one 
hand, the constraint error rate is less than when no constraint correction is employed, but 
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is not held constant as in the previous method. 

On the other hand, the attitude error rate (as 
defined previously) is not constant, but in- 
creases so that the results must be presented 
in terms of rate of change of attitude error 
rate per cycle, This higher order increase in 
attitude error indicates that normalization is 
not an effective technique for improving quater- 
nion kinematic simulations. The results of this 
study are given in Figure 7. 

DERIVATIVES’ CONSTRAINT 

At this point in the study, all the techniques for reducing the constraint error discussed in the 
literature have been studied. An attempt was made to develop a new technique which may serve 
the same purpose and be practicable also. 

The most natural idea seems to be to select a function f and to approximate the equations for e 
by the equations representing (e + f ) for all e T s. This function f may serve as a feedback term to 
force the constraint error toward zero; however, it must not disturb the stability of the system. 

If each equation in Equation 16 is squared and the resulting equations added together, simpli- 
fying the results gives 


CONSTRAINT ERROR RATE (eye' 1 ) 

RATE OF CHANGE OF ATTITUDE ERROR RATE 

( eye/ eye/ eye ) 



Figure 7— Quaternions with normalization. 


e Q 2 + 42 + e 2 2 + 4 3 2 


( e o 2 + e i 2 


+ e 2 + e 2 
2 3 



(29) 


where 


+ CO 2 


If no constraint error enters the solution, Equation 29 is simply 



(30) 


thus a new constraint equation has been derived based on the derivatives of the quaternions. We 
can define an error e d based on this constraint equation: 


2 + 42 + e 2 + 


co 2 
M 


(i-O T 


co 2 
M 

4 


(31) 
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where 


e 


d 



(32) 


It is seen that the error e d is caused by the constraint error e, which affects the derivatives e , 
e l , e 2? and e 3 directly. The function f may be chosen to take the form ee 0 , ee J? ee 2 , and ee 3 for 
each one of the equations in Equation 16, respectively. The set of differential equations thus as- 
sumes the form 


(*o)t = + 

(*l)« = T ( e 0^M 1 “ e 3 & M2 +e 2 a ' M3 ) t + “( e >. (®l)t-h ’ 
( 4 2 )t = i ( e 3 W M 1 +e 0 C "M 2 " e i £ ' , M 3 ) t + M ( £ )t ( 4 2 )t-h • 






M(«)t ( 4 3 ) 


t- h ’ 


(33) 


where M is a constant. The subscript used outside the parenthesis indicates the time at which the 
quantity inside is being evaluated. In dealing with these equations, the constraint error, the qua- 
ternion components, and their corresponding derivatives are computed once at the end of every 
time interval. Their substitution into Equation 33 prepares the derivatives for the next interval. 

The procedure for finding the stable value of M is similar to that used for K as described in 
the algebraic constraint. When & and h are not normalized, then the maximum stable value of M 
associated with the minimum constraint error as obtained from computer solutions is, not only a 
function of step size, but also a function of angular speed. It was demonstrated by computer solu- 
tions that the a> u and h associated with optimum M satisfy a reciprocal relationship; that is, the re- 
sult is the same if h is multiplied by a constant, and is divided by the same constant. This is 
demonstrated mathematically as follows: 


a 



L = 1, 2, 3 , 


and 


r = Bt , 
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where A and B are two constants. Then 


„ de j _ d _!i dr - 

e j dt ~ dr dt Be. , 

j 3 0, 1, 2, 3 - 


Hence, 

e 0 = 'H ( e i n M 1 +e 2°M 2 +e 3 n M 3 ) + M£e 0' • 

*1 = II ( e 0 n Ml “ e 3 n M 2 +e 2 n M 3 ) + M£e i' ’ 

e 2 = ll" ( C 3 Q M 1 +e 0 Q M 2 ‘ C i n M 3 ) + M£e 2' > 

e 3 ‘ - II ( e 2 U M, - e !°M 2 " e 0 U M 3 ) + Mte 3' ' (34) 


If we set A = b, then the equations have the same form as the original set of equations (Equa- 
tion 33); they will be satisfied by the same solution and thus have the same optimum M. However, 
for the original set, the step size is 


At - h , 


and for the second set 


A'/ BAt H , 


or 


M • I * (35) 

Thus H ■ a is equal to h * and the product h • a, can be treated as a single parameter. Next an ex- 
pression was sought relating M and h ■ w in a simple form. As before, the practical interest centers 
mainly on the values of M such that the constraint error is diminished to the limit for which the 
system remains stable. 
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Extensive review of computer solutions suggested the expression 


M • ( h * j 2 < 2 /tt 2 


for (h • w M ) in cycles, or 

M • ( h -"m) 2 < 8 (36) 

for (h * o> M ) in radians. 

As previously, an attempt is made to explain Equation 36 analytically. Note first that an inde- 
pendent stability condition | M( e ) 1 1 < 1 can be established by inspection (e.g., for w = 0) from Equa- 
tion 33. To establish Equation 36, the first equation of Equation 33 is multiplied by (e ) t , the sec- 
ond by (e 1 ) t , the third by (e 2 ) t , and the fourth by (e 3 ) t . When these four equations are summed 
and simplified, 


( e 

\ o 0 


+ e e, 4 e„ + e e 0 
1 1 2 2 3 3 


)t = M(e) t [( e o)t( 4 o)t-h + ( e i)t( 4 ,)t-h + ( e 2 )t( 4 2 )t-h + ( e 3 )t(s)t-h] • (37) 


That is, 



M ( e )t [( e o)t(®o)t-h + ( e i)t(^)t-h + ( e 2 )t (e 2 )t-h + ( e 3 )t (® 3 )t-h] • 


(38) 


The purpose of this technique is to constrain the derivatives such that the difference of &.• 2 /4 and 

M 

the sum of the square of the quaternion components' derivatives decrease. For simplicity, (e) t 
may be written approximately as 


( e > t = ( e )t- h + h ( 4 ) t -h 


(39) 


Substituting Equation 39 into Equation 38 gives 


(AO t 

2h~ = M ( £ >t 


< Ae )«- h W , ' 

2h 4 U ^ e )t-h) 


( 40 ) 


That is, 


M(O t hW 

2 (l-( £ )t- h ) = M(e) t (A€) t . h - (A6) t 
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Because (e) t r 0, and (e) t . h « 1, we may write 


Mh 2 ^ 2 (A€) t _ h (Ae) t 

2 * M ( £ )t (€) t (S) t 


( A O t - h 


1> 

M(O t (e) - 

+ 

( £ ) t 


Since M(e) t < 1 for stability, 


M h 2 co* 

( A O t -h 


( A O t 

2 < 

(Ot 

+ 

(Ot 


Near an oscillatory stability boundary, 


( a O t _ h - " (Ae) t % (2 £ ) t 


giving 


M h 2 w 2 

M 


< 4 


or 


(41) 

(42) 


Mh 2 ^ <8 • 

To show that this condition includes the entire stable region, consider the value of M such that 
Mh 2 oj 2 > 8 and assume that at time t -h the constraint error is small enough so that I M( e ) t _ h < 1. 

M 1 1 

By using Equation 42, we have 


<TO t - h 


( A O t 

Mh 2 o> M 2 

M(O t (e)t 

+ 

(Of 

“ 2 


> 4 . 


This implies that either 


(Ae) t 

(O t 


> 2 


or 


< A O t - h 

M < e >t tot 


> 2 ■ 
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STEP SIZE (cycle) 


The first of these expressions implies that after time t 
fore unstable. 


and the solution is there- 


, (A6) t >2 (e) t , 


If the second inequality holds, we can write 


tr-h 

M( € ) t ^ ^ “ M( £ ) t 




(O t 


= M(e )t - M(6) t _ h 


This indicates that 


M(e) t -M(e) t _ h > 2 


or 


|M(O t | + |M(O t _ h | > 2 ■ 

However, it was assumed originally that 


|M(Ot- h | < 1 ’ 

thus 


M(e) t | > 1 


CONSTRAINT ERROR 

ATTITUDE ERROR RATE (cycle/cycle) 


-8 -7 -6 -5 -4 

10 10 10 10 10 



0 20 40 60 80 100 120 140 160 180 200 

VALUE OF M ABOVE WHICH TECHNIQUE IS UNSTABLE 

(<*>= In RAD) 

Figure 8— Quaternions with derivative constraint. 


driving the solution unstable. 

Since M depends on ( ha> M ) 2 , it is natural to 
ask what is the stable value of m if h is held 
constant, but a> M is allowed to vary 0 In this 
case, must be replaced by its maximum for 
the same h in Equation 36 to obtain the value of M 0 

Figure 8 represents the results of this 
technique. By comparison with the results for 
the algebraic constraint shown in Figure 6, it 
can be seen that both the attitude error rate 
and the constraint error are reduced. 


20 




■ II I 


CONCLUSIONS 

This study has produced several significant results. The findings developed in the study of the 
benchmark problem are rather general and have application to all rotational kinematics simula- 
tions. The conclusions of this study are summarized as follows. 

1. Euler angles have visual appeal, but attitude error increases as the singularity is approached. 

2. The quaternion method is superior for accuracy and speed. 

3. The fixed-step integration technique appears superior to the variable-step technique in 
simplicity and speed for simulations involving reasonable frequency variations. The relative merits 
for simulations involving solutions with widely varying frequency and/ or amplitude can be inferred 
from Figures 3 and 4, but are strongly problem-related and cannot be generalized. 

4. The error analyses indicate how much error may be expected in the solution for a pre- 
assigned step size. The results presented graphically will also assist in estimating simulation 
speeds. 

5. In the quaternion method, the algebraic constraint technique is useful for reducing the con- 
straint error, but it has no effect on the attitude error rate. 

6. The newly proposed derivatives’ constraint technique is not only as effective as the alge- 
braic constraint technique in constraint error correction, but also diminishes the attitude error 
rate to a small constant value. A comparison can be made by a close examination of Figures 6 ‘and 8. 

7. The normalized constraint technique provides a small reduction in constraint error buildup, 
but it fails to stop the increasing of the attitude error rate. Consequently, this technique is not 
efficient and has less value in application. 

8. The numerical stability condition for each technique has been discussed and determined. 

When algebraic constraint is employed, the stability condition is h • K < 1, where h is the integra- 
tion step size, and K is a positive constant used in implementing the constraint equation. When 
derivatives’ constraint is used, the stability condition is Mh 2 u? < 8, where is angular speed (rad/ 
sec), and M is a positive constant used in implementing the constraint equation. 
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